#WIC co-op eWIC new diff in diff pictures
#last modified: 6 May 2025
#last modified by: Charlotte Ambrozek
#this do file takes estimates from the csdid2 package for groups, years, and subsamples with chain by state timing and constructs graphs as needed
library(haven)
library(tidyverse)
library(RColorBrewer)
library(ggpmisc)

event <- read_dta("./analysis/output/csdid_tip_chainstatetime_event.dta") 

picture <- filter(event, state == "all") %>%
  group_by(outcome, subsample, controls, ev_yr) %>%
  mutate(subsample = gsub("all", "All", subsample)) %>%
  mutate(subsample = gsub("chain", "Chain", subsample)) %>%
  mutate(subsample = gsub("indep", "Independent", subsample)) %>%
  mutate(controls = gsub("no", "No", controls)) %>%
  mutate(controls = gsub("yes", "Yes", controls))

# make three pictures:
#one that just shows all
#one that compares chain and indep
#one that compares all to highwicelig

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All" & controls == "No")
allpre <- filter(picture, subsample == "All" & controls == "No", ev_yr < 0)
allpost <- filter(picture, subsample == "All" & controls == "No", ev_yr >= 0)

all_plot <- ggplot(all, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allpre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allpost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allpost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year") +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_chainstatetime_auth_all.png", plot = all_plot, width = 6.5, height = 4, units = "in")

allcontrols <- filter(picture, subsample == "All" & controls == "Yes")
allcontrolspre <- filter(picture, subsample == "All" & controls == "Yes", ev_yr < 0)
allcontrolspost <- filter(picture, subsample == "All" & controls == "Yes", ev_yr >= 0)

allcontrols_plot <- ggplot(allcontrols, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = allcontrolspre, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allcontrolspre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = allcontrolspost, aes(x = ev_yr, y = avgb), colour = "gray70") + 
    geom_ribbon(data = allcontrolspost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate",
    x = "Event Year") +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_chainstatetime_auth_all_controls.png", plot = allcontrols_plot, width = 6.5, height = 4, units = "in")

chain <- filter(picture, (subsample == "Chain") & controls == "No") %>%
    group_by(outcome, ev_yr) 
chainpre <- filter(chain, ev_yr < 0) %>%
    group_by(outcome, ev_yr) 
chainpost <- filter(chain, ev_yr >= 0) %>%
    group_by(outcome, ev_yr) 

# Create the chain plot
chain_plot <- ggplot(data = chain, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = chainpre, aes(y = avgb, x = ev_yr), color = "gray70") + 
    geom_ribbon(data = chainpre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = chainpost, aes(y = avgb, x = ev_yr), color = "gray70") + 
    geom_ribbon(data = chainpost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate", x = "Event Year") +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0), legend.position="none")
  
ggsave("csdid_tip_chainstatetime_auth_chain.png", plot = chain_plot, width = 6.5, height = 4, units = "in")

# make chain plot with controls
chaincontrols <- filter(picture, (subsample == "Chain") & controls == "Yes") %>%
    group_by(outcome, ev_yr) 
chaincontrolspre <- filter(chaincontrols, ev_yr < 0) %>%
    group_by(outcome, ev_yr) 
chaincontrolspost <- filter(chaincontrols, ev_yr >= 0) %>%
    group_by(outcome, ev_yr) 

# Create the chain plot
chaincontrols_plot <- ggplot(data = chaincontrols, aes(x = ev_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_vline(xintercept = 0, lty = 3, colour = "red") + 
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = chaincontrolspre, aes(y = avgb, x = ev_yr), color = "gray70") + 
    geom_ribbon(data = chaincontrolspre, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    geom_line(data = chaincontrolspost, aes(y = avgb, x = ev_yr), color = "gray70") + 
    geom_ribbon(data = chaincontrolspost, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate", x = "Event Year") +
    theme_light() +
    scale_x_continuous(breaks = seq(-4, 4, 2)) +
    theme(plot.caption = element_text(hjust = 0), legend.position="none")
  
ggsave("csdid_tip_chainstatetime_auth_chaincontrols.png", plot = chaincontrols_plot, width = 6.5, height = 4, units = "in")

# make processor graphs
# load in implementation years to use when figuring out implementation order
ewic <- read_dta("./data/cleaned/ewic_rollout.dta") %>%
  select(c(state, ev_year, ct_fips)) %>%
  distinct()

tip <- read_dta("./data/cleaned/tip_auth_sq.dta") %>%
    select(c(state, ct_fips)) %>%
    distinct() %>%
  rename(st = state)

process <- left_join(tip, ewic, by = "ct_fips") %>%
  filter(!((st == "KY" & state == "Tennessee") | (st == "NH" & state == "Maine"))) %>%
    mutate(processor = case_when(st %in% c("VA", "CT", "NH", "NY", "MI", "OK", "RI", "AL", "GA", "MS", "SC", "TN", "IN") ~ "Conduent",
                            st %in% c("NM", "TX", "ME", "MD", "NJ", "PA", "NC", "IL", "AR", "LA", "UT", "MT") ~ "Solutran", #note that solutran implemented wyoming in 2002 before new mexico or texas
                            st %in% c("MA", "WV", "FL", "KY", "WI", "OR", "DE", "DC", "IA", "MI", "MN", "AZ", "CO", "KS", "MO", "NE") ~ "FIS CDP")) %>%
    group_by(st) %>%
    mutate(ev_year = min(ev_year)) %>%
    mutate(avg_ev_year = mean(ev_year)) %>%
    select(c(st, ev_year, processor, avg_ev_year)) %>%
    rename(state = st) %>%
    distinct()

simple <- read_dta("./analysis/output/csdid_tip_chainstatetime_simple.dta") 

processor_all <- filter(simple, state != "all" & subsample == "all" & outcome == "auth" & controls == "no") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
    filter(!is.na(processor) & !is.na(se))

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")

# Define the shapes for each chain status
shapes <- c(16, 17, 18)

# Add new columns for color and shape mapping
processor_all$color <- factor(processor_all$processor)
processor_all$shape <- factor(processor_all$processor)
# Create the chain plot
processor_all_plot <- ggplot(data = processor_all, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate", x = "Processor order") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 2)

ggsave("csdid_tip_chainstatetime_auth_processor.png", plot = processor_all_plot, width = 6.5, height = 4, units = "in")

processorcontrols_all <- filter(simple, state != "all" & subsample == "all" & outcome == "auth" & controls == "yes") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
    filter(!is.na(processor) & !is.na(se))

# Define the colors for each chain status
colors <- brewer.pal(3, "Set2")

# Define the shapes for each chain status
shapes <- c(16, 17, 18)

# Add new columns for color and shape mapping
processorcontrols_all$color <- factor(processorcontrols_all$processor)
processorcontrols_all$shape <- factor(processorcontrols_all$processor)
# Create the chain plot
processorcontrols_all_plot <- ggplot(data = processorcontrols_all, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate", x = "Processor order") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 2)

ggsave("csdid_tip_chainstatetime_auth_processor_controls.png", plot = processorcontrols_all_plot, width = 6.5, height = 4, units = "in")

processor_chain <-  filter(simple, state != "all" & subsample == "chain" & outcome == "auth" & controls == "no") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_chain$color <- factor(processor_chain$processor)
processor_chain$shape <- factor(processor_chain$processor)
# Create the chain plot
processor_chain_plot <- ggplot(data = processor_chain, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, chain stores", x = "Processor order") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_chainstatetime_auth_processor_chain.png", plot = processor_chain_plot, width = 6.5, units = "in")

processor_chain_controls <-  filter(simple, state != "all" & subsample == "chain" & outcome == "auth" & controls == "yes") %>%
    left_join(process, by = "state") %>%
    group_by(processor) %>%
    arrange(processor, ev_year, avg_ev_year) %>%
    mutate(n = row_number()) %>%
  filter(!is.na(processor) & !is.na(se))

# Add new columns for color and shape mapping
processor_chain_controls$color <- factor(processor_chain_controls$processor)
processor_chain_controls$shape <- factor(processor_chain_controls$processor)
# Create the chain plot
processor_chain_controls_plot <- ggplot(data = processor_chain_controls, aes(x = n, y = b, color = color, shape = shape)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_smooth(method = lm, se = FALSE) +
    stat_poly_eq(formula = y ~ x,
                output.type = "numeric",
                mapping = aes(label = sprintf("Slope: %s",
                            c(formatC(round(stat(coef.ls)[[1]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[2]][[2, "Estimate"]], digits = 3)),
                              formatC(round(stat(coef.ls)[[3]][[2, "Estimate"]], digits = 3))))),
                label.x = 3.5,
                label.y = -0.05) +
    geom_hline(yintercept = 0, lty = 2) +
    scale_color_manual(values = colors, name = "Processor") +
    scale_shape_manual(values = shapes, name = "Processor") +
    labs(y = "ATT estimate, chain stores", x = "Processor order") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    facet_wrap(~ processor, ncol = 1)

ggsave("csdid_tip_chainstatetime_auth_processor_chain_controls.png", plot = processor_chain_controls_plot, width = 6.5, units = "in")

### make treatment group specific graphs
group <- read_dta("./analysis/output/csdid_tip_chainstatetime_group.dta") 

picture <- filter(group, state == "all") %>%
  group_by(outcome, subsample, gr_yr, controls) %>%
  mutate(subsample = gsub("all", "All", subsample)) %>%
  mutate(subsample = gsub("chain", "Chain", subsample)) %>%
  mutate(subsample = gsub("indep", "Independent", subsample)) 

avg <- filter(picture, gr_yr == "Average") %>%
    ungroup() %>%
    select(!c(gr_yr)) %>%
    mutate(avgb = b, avgll = ll, avgul = ul) %>%
    select(c(outcome, subsample, state, controls, avgb, avgll, avgul))

picture <- filter(picture, gr_yr != "Average") %>%
  mutate(gr_yr = as.numeric(gr_yr))
picture <- left_join(picture, avg, by = c("outcome", "subsample", "state", "controls"))

rm(avg)

# make three pictures:
#one that just shows all
#one that compares chain and indep
#one that compares chain and indep with different timings

setwd("./analysis/output/graphs")

all <- filter(picture, subsample == "All" & controls == "no")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented") +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_chainstatetime_auth_all_group.png", plot = group_plot, width = 6.5, height = 4, units = "in")

chain <- filter(picture, controls == "no" & (subsample == "Chain"))

# Create the chain plot
chain_plot <- ggplot(data = chain, aes(x = gr_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = chain, aes(y = avgb, x = gr_yr), color = "gray70") + 
    geom_ribbon(data = chain, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2))
  
ggsave("csdid_tip_chainstatetime_auth_chain_group.png", plot = chain_plot, width = 6.5, height = 4, units = "in")

all <- filter(picture, subsample == "All" & controls == "yes")

group_plot <- ggplot(all, aes(x = gr_yr, y = b)) +
    geom_point(size = 3) +
    geom_ribbon(aes(ymin = all$avgll[1], ymax = all$avgul[1]),
              alpha = 0.25,
              color = "gray90") +
    geom_linerange(aes(ymin = ll, ymax = ul),
                alpha = 0.7,
                size = 0.8) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_hline(yintercept = all$avgb[1], color = "gray70") + 
    labs(y = "ATT estimate",
    x = "FY WIC EBT First Implemented") +
    theme_light() + 
    scale_x_continuous(breaks = seq(2006, 2016, 2)) +
    theme(plot.caption = element_text(hjust = 0))
  
ggsave("csdid_tip_chainstatetime_auth_all_group_controls.png", plot = group_plot, width = 6.5, height = 4, units = "in")

chain <- filter(picture, controls == "yes" & (subsample == "Chain"))

# Create the chain plot
chain_plot <- ggplot(data = chain, aes(x = gr_yr, y = b)) +
    geom_linerange(aes(ymin = ll, ymax = ul),
                   alpha = 0.7,
                   size = 0.8) +
    geom_point(size = 3) +
    geom_hline(yintercept = 0, lty = 2) +
    geom_line(data = chain, aes(y = avgb, x = gr_yr), color = "gray70") + 
    geom_ribbon(data = chain, aes(ymin = avgll, ymax = avgul),
                   alpha = 0.25,
                   colour = "gray90") +
    labs(y = "ATT estimate", x = "FY WIC EBT First Implemented") +
    theme_light() +
    theme(plot.caption = element_text(hjust = 0), legend.position="none") + 
    scale_x_continuous(breaks = seq(2006, 2016, 2))
  
ggsave("csdid_tip_chainstatetime_auth_chain_group_controls.png", plot = chain_plot, width = 6.5, height = 4, units = "in")